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Abstract 

Intensity distributions of linear wave fields are, in the high frequency limit, often approximated 
in terms of flow or transport equations in phase space. Common techniques for solving the flow 
equations for both time dependent and stationary problems are ray tracing or level set methods. 
In the context of predicting the vibro-acoustic response of complex engineering structures, reduced 
ray tracing methods such as Statistical Energy Analysis or variants thereof have found widespread 
applications. Starting directly from the stationary Liouville equation, we develop a boundary el- 
ement method for solving the transport equations for complex multi-component structures. The 
method, which is an improved version of the Dynamical Energy Analysis technique introduced re- 
cently by the authors, interpolates between standard statistical energy analysis and full ray tracing, 
containing both of these methods as limiting cases. We demonstrate that the method can be used 
to efficiently deal with complex large scale problems giving good approximations of the energy 
distribution when compared to exact solutions of the underlying wave equation. 

Keywords: Statistical energy analysis, high-frequency asymptotics, Liouville equation, boundary 
element method 



1. Introduction 

Many phenomena in physics and engineering can be accurately described or modelled in terms 
of linear wave equations. A large range of numerical methods has been developed for solving 
wave problems, often making heavy use of the linearity of the underlying PDEs leading in general 
to large, finite matrix equations. Popular tools include finite element or finite volume methods, 
boundary element methods and various spectral methods. There are, however, basic limitations 
when approximating the solutions of the wave equations directly: the size of the associated linear 
system increases with decreasing wavelength and numerical schemes become inefficient when the 
local wavelengths are orders of magnitude smaller than typical dimensions of the physical system. 

In order to overcome this limitation, a range of high-frequency/ short wave length methods 
have been devised such as the WKB or eikonal approach or semiclassical methods developed in the 
context of quantum chaos, see [TJ [H |3] for overviews. In particular, phase information is obtained 
by solving the Hamilton- Jakobi equations for the action fields, whereas amplitude information 
is calculated using transport equations. The actual methods can be quite intricate demanding 
an intimate knowledge of the classical phase space dynamics given by the Hamilton equations of 
motion for the underlying ray or classical dynamics. If one is prepared to omit phase information 
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altogether, the approximations amount to solving the transport or Liouville equation for phase 
space densities. 

In this paper, we will focus on efficient numerical methods for solving the Liouville equation 
(LE) for stationary problems. We will show how the solutions of the LE are related to both ray 
tracing (see for example, [HOIS]) and statistical energy analysis (SEA) [7J [HI IH], a popular method 
in the engineering community for solving vibro- acoustic problems in the high-frequency limit. 

The Liouville equation is an example of a PDE describing a conservation law |10| . Typical 
solution methods are based on characteristics or ray tracing. Ray tracing is the method of choice 
when considering propagation over relatively short times or length intervals. Prime examples include 
weather forecasting, where data updates lead to regular re-initialisations of the density distributions 
[TT] . or room acoustics, where only very few reflections need to be considered [5]. More sophisticated 
methods related to tracking the time-dynamics of interfaces in phase space such as moment methods 
and level set methods [121 113 IH] have been developed only relatively recently, finding applications 
in acoustics, seismology and computer imaging. For an excellent recent overview, see Runborg [2]. 
Numerical methods based on solving the time dependent Liouville equation directly on a mesh using 
finite volume or transfer operator methods, also referred to as Ulam methods, have been considered 
in [T5] and more recently in [THl [UJ . Ray based and tracking methods often become inefficient when 
considering frequency domain wave problems in bounded domains, for example, determining the 
wave field in a finite size cavity driven by a continuous, monochromatic excitation. Here, multiple 
reflections of the rays and complicated folding patterns of the associated level-surfaces often lead 
to an exponential increase in the number of branches that need to be considered. 

In the engineering literature, these problems have been circumvented by subdividing the struc- 
ture into a set of substructures. Assuming ergodicity of the underlying ray dynamics and quasi- 
cquilibrium conditions in each of the subsystems, that is, assuming that the density in each subsys- 
tem is approximately constant, greatly simplified SEA equations based only on coupling constants 
between subsystems can be set up. The disadvantage of SEA is that the underlying assumptions 
are often hard to verify a priori or are only justified when an additional averaging over 'equivalent' 
subsystems is considered. These shortcomings have been addressed by Langley [TH1 EH] and more 
recently in a series of papers by Le Bot [213 [211 122] ■ A computational tool based on a Frobenius- 
Perron operator approach systematically interpolating between SEA and full ray tracing has been 
introduced first in |23j and further analysed in |24| . Its name, Dynamical Energy Analysis (DEA), 
points at the similarities with SEA but stresses at the same time the importance of dynamical cor- 
relations present in the ray-densities. Thus, in DEA, the solution in a subsystem of the structure is 
typically represented by a vector; this is in contrast to SEA where the solution in each subsystem 
is represented by a single number, the mean ray density in each subsystem. 

The implementation of DEA as presented in [23J [21] corresponds to a spectral boundary integral 
method; the integral equations are expanded in orthogonal basis sets and the resulting matrix 
equations are solved for the coefficients of the basis expansion of the solution. It has been shown 
in |24j that the method works even in situations where a naive application of SEA fails. However, 
spectral methods can lead to slow convergence for the matrix elements in the DEA matrix, in 
particular if the boundary is non-smooth. In this paper, we show that using a boundary element 
method (BEM) for the spatial variable and a basis function expansion in the momentum component 
leads to large efficiency gains. We will furthermore present a more in depth derivation of the actual 
boundary integral formulation. 

The paper is structured as follows: in Sec. [2j we sketch the derivation of the Liouville transport 
equation from the underlying wave equation; we will then derive the boundary integral formulation 
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in Sec. [3] and give a description of the numerical implementation of DEA using the BEM/ basis 
function approach in Sec. |4j Numerical results will then be presented in the Sec. [5j 



2. From wave dynamics to phase space densities 

2.1. A time dependent formulation 

In the following, we sketch the connection between linear wave equations and the Liouville 
equation via the Eikonal ansatz. We develop the theory in the time domain first and move to a 
stationary formulation in the next section. We will restrict the discussion to the Helmholtz equation 
of acoustics; other wave problems such as electromagnetic waves, linear elasticity or Schrodinger's 
equation can be treated along the same line with the obvious modifications. For a more in depth 
discussion of the Eikonal approach, see for example [5j ■ A nice derivation of the connection between 
wave transport equations and the Liouville equation based on Wigner transformation techniques 
has been given in [2"5] . 

We start from the wave equation for linear acoustics, that is, 



seeking solutions u(r, t) in some domain f2 with boundary T (u may express acoustic pressure in 
appropriate units). Furthermore, c = c(r) > is the local wave speed. We assume in what 
follows that the wave equation is not explicitly time dependent and that the walls act as passive 
elements, that is, we have Dirichlet or Neumann boundary conditions, u(s) — or du/dn(s) = 
for s € r. Solutions are then determined uniquely in terms of initial conditions u(r, 0) = uo(r) and 
d t u(r,0) — iti (r) at t = 0. 

We now make the usual Eikonal ansatz 



with real valued functions A — A(r, t), <j> = (f>(r, t) and ui acting as the large parameter here. Plugging 
(pi) into M and collecting terms of the order uj 2 , we arrive at the Eikonal equation 




(1) 



u = Ae 



(2) 




(3) 



Terms of order u> yield a transport equation 




- 2c 2 \7A 



- c 2 AA<f> - 2cV • (cA\7(f>) = 0. 



(4) 




(5) 
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with "velocity field" 

v (r,i) = - c 2 V0/^ (6) 

determined through the phase function <fi(r, t) alone. 

Solutions of the Eikonal equation (J3| can readily be found using characteristics. Without re- 
stricting the discussion we choose the "+" sign in (J3j) ; see also comments further down in the text. 
Setting 

H(r,p)=c\p\ with momentum p — V(f>, (7) 

one obtains solutions 4>(r,t) along trajectories in phase space X(t) — (r(t),p(t)) by solving Hamil- 
ton's equations of motion, that is, 



V p H = (8) 



P_ 
>l 

-V r i? = |p|Vc 



Note that in deriving the Eikonal Eq. ^ and transport Eq. (J4J) we omitted terms of order O(ui ); 
neglecting the coupling terms between Eqs. ([3]), Q leads to first order differential equations which 
can formally be solved by using only one of the initial conditions, that is, either uo(r) or u\(r). 
The ODEs Q are thus solved by using, for example, initial conditions X(0) — (ro,po == ^4>o{. r o))- 
A solution for the phase </> is then obtained by integration, that is, 

<K?",*) = <t>o(ro) + / Prdr 
Jo 

where the integral is taken along X(t) = ^(Xo). Here ^(Xo) is the phase space flow map obtained 
from the pair of ODEs (J8|. In particular we integrate over trajectories starting at t = at some 
point r with momentum p — V4>o(r ) and reaching the final point r at time t. Here, 4>o is t ne 
initial phase corresponding to the initial conditions uo( r o) = ^o(^o) exp(iw0o( r 'o)) at t = 0. In 
general, there will be several such initial points from which a trajectory reaches the point r at time 
t. That is, the solution to the eikonal equation |3| at a point r has a multi-sheet structure with a 
set of N — N(r, t) branches, each giving rise to a phase 

4>j{r, t) = M^j) + f PjTjdr, j = 1, N . (9) 
Jo 

The integrals are taken over paths Xj, j = 1, N, reaching r in time t with Xj(0) = (rj, V0o( r j))- 
It is this multi-sheet structure which makes solving Eq. ^ difficult in practice [2]. In particular, 
the number of branches N(t) can increase rapidly with t and for chaotic ray dynamics N(t) actually 
increases exponentially pQ. 

We would like to remark that choosing the "— " sign in ^ leads to a sign-change in ^ and 
gives time-reversed ray trajectories obtained through the time transformation t — > —t. The complete 
solution to the Eikonal equation is thus given by the set of characteristics travelling both in positive 
and negative time. 

We turn now to the solution of the transport Eq. ([5]); this equation is driven by the solution 
4>{r, t) of the Eikonal Eq. ([3]). It thus suffers in principle from the same multi-valuedness problem as 
the original equation. Note, however, that the solutions of the ray equations (pi in phase space are 
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determined uniquely for a given initial condition; it is the projection onto r - space which leads to 
the multi-valuedness of the phase solution. It is thus advantageous to solve the transport equation 
in phase space. The corresponding equation is the Liouville equation (LE) 

f t +{H,p} = (10) 

where {• , •} are the Poison brackets, and H(r,p) is the classical Hamilton function Q. The phase 
space density 

N(t) 

p(r,p,t) = ^p J -(r,i)£(p-V^(r,i)) (11) 



indeed solves Eq. ( |l0[ ) if the branches 4>j ( r , t) and (jj (r, t) obey the Eikonal and transport equations 
([3]) , (|5J(, respectively. (This result can be obtained by inserting ( 11 ) into ( 10 1, see also [2J.) A formal 



solution of the LE is obtained in terms of a Frobenius-Perron operator 

p{X,t)=C t [ P ](X)= f S(X-^ t (Y))p(Y,0)dY = p^- t (X),0), (12) 



where the last relation is valid for Hamiltonian flows; this relation is used extensively in the ray 
tracing approach. 

The solution of the wave equation is now formally obtained by summing over the solution 
branches in ([9]), that is, 

N 

u(r,t) = y £ j \A j {r,t)\e iu<t> ^- iv ^ (13) 
3=1 

where Vj are Maslov indices [T] due to caustics along the ray Xj(t). As we are not interested in 
phases in what follows, we will not discuss this issue any further here. An asymptotic treatment 
of the wave equation using the Eikonal ansatz needs to keep track of the link between phases and 
amplitudes along the different solution branches. Ray tracing methods are often the only way 
forward when finding solutions using the Eikonal ansatz. This may be an efficient method for short 
times and in cases where only a few reflections occur. However, the rapid increase in the number 
of paths makes ray tracing methods quite cumbersome for large times. The advantages behind 
the original idea of the Eikonal ansatz, namely replacing a wave equation with highly oscillatory 
solutions by differential equations for slowly varying quantities such as the amplitude and phase, 
are lost. 

In many applications it is, however, sufficient to find the solution of the transport equations for 
a given set of initial conditions. The density p(X, t) contains a great deal of information about the 
wave solution, and this may often be sufficient for the required purpose. Classical ray tracing, that 
is ray tracing without tracking phase information, is indeed used extensively in room acoustics. 

2. 2. Stationary formulation 

In what follows, we will consider wave problems with time-harmonic forcing at a frequency ui. 
We can then separate out time by setting 

u(r,t) = G(r)e tut . 
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Let us consider problems with a forced excitation by a point source. The corresponding inhomoge- 
neous wave equation then has the form 

(V ■ (c 2 V) + w 2 ) G(r) = -c 2 S{r - r ) , (14) 

with the source point at r = r . The solution G(r, r^) is the Green function and arbitrary force 
distributions can be recovered. 

For an Eikonal approach, we proceed by splitting the solution into a homogeneous and an 
inhomogeneous (or direct) solution in the form 

G(r,r ) = G (r,r ) +G h (r,r ), 

where Go is the free Green function and Gh solves the homogeneous Helmholtz Eq. 

(V-(c 2 V)+u 2 )G h {r)=0 (15) 

with appropriate boundary conditions. (For Dirichlet BC, that is G(s) = 0, we obtain Gh{s) — 
—Gq(s) for s € r). One finds for two-dimensional problems with c = const, 

G (r,r ) = ~HZ(u>\r-r \/c), (16) 

where Hq is the 0th order Hankel function of the first kind. Proceeding with an Eikonal ansatz, we 
set Gh = Aexp(iujip)] we have then 

dA 

cb(r,t) = ib(r)+t and — =0. (17) 
at 

The Eikonal and transport Eqs. pH), ([5| reduce to 

c|VV>| = 1 (18) 

and 

V • (vp) = with v(r) = -c 2 Vt/>. (19) 
We write the boundary conditions in the form 

G (s) = A (s)e^ o(s) forseT. 



In order to obtain an approximate solution of ( |18[ ) at a point r £ £1 one needs to find those rays 
starting from a boundary point s £ V with initial conditions Xq = (s,p = VV'o( s )) which pass 
through the point r following the equations of motion The phase function ijj is then obtained 
by integrating along this trajectory with an initial condition ip(s) — ipo(s); the phase function ip 
again has a multi-sheet structure and there are in general infinitely many sheets intersecting with 
the plane r = const in phase space. (A phase space description of the sheets can be given in terms 
of Lagrangian manifolds, see for example [26].) The various branches are calculated according to 



V'iO") = Msj) + pdr 
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integrating along a trajectory with initial condition Xj = (sj,pj = V - 0o(sj))- The transport 
equation ( 19 1 is solved with initial conditions 



j (s) = -<?A\a) 



, dip(s) 
dn 



(20) 



where jo{s) is the incoming flux normal to the boundary. The total solution is then approximativcly 
obtained as 

G(r) = GoO) + l^-Mle^M-^/ 2 (21) 

3=1 

where pj are again Maslov phases and \Aj \ is the amplitude transported along the trajectory from 
Xj to r; for alternative derivations of Eqs. (13 1, (21 1, see pQ. 

In applications, one is often interested in wave intensities (such as energy densities); the contri- 
butions arising from reflections from the boundary have the form 

I(r) = |G,(r)| 2 = ]T \A,A r \e^-^ = £>(r) + £ \A 3 A r \e^~^ . (22) 

Neglecting the oscillatory terms, which is justified after averaging over a small U) range or over an 
ensemble of similar systems, one arrives at the ray tracing approximation 



7(r) « p( r ) = ^Pj(r) = J P( r iP) d P 



(23) 



where p(r,p) is a solution of the stationary Liouville equation with 

{H,p} = 



(24) 



and boundary conditions given by ( 20 ) together with Snell's law of reflection. Note that the total 
intensity |G(r) | 2 for r 6 Q may be estimated simply by adding |Go(?')| 2 to p{r). 

It is clear from the description above that a ray tracing approach is even more cumbersome 
in the stationary case than in the time-dependent case. The stationary solution is constructed 
from infinitely many branches which originate from reflections at boundaries. If we neglect phase 
information, however, we can reconstruct the phase space densities using the Frobenius-Perron 
operator, Eq. (12 1 [231 IMj- In the next section we will introduce an efficient algorithm for deter- 
mining solutions of the stationary Liouville equation even in the presence of multiple reflections, 
thus overcoming the limitations of the ray tracing approaches sketched above. 



3. Boundary integral formulation of the Liouville equation 



In order to solve the stationary Liouville equation ( 24 1 , the ray tracing integral formula, Eq. 
( 12 ), is rewritten in a boundary integral form using a boundary mapping technique. The boundary 
mapping procedure involves first mapping the ray density emanating continuously from the source 
onto the boundary T and then constructing a Frobenius-Perron operator for the boundary map. 

We now give an explicit derivation of the initial density on the boundary arising from Eq. ( 14 ) 
for Dirichlet boundary conditions. The ray density emanating from a source point r G is given 



7 



by the square of the amplitude of the free space Green function. For C R 2 , it is given by (16) 
and we have the following high frequency asymptotic form 

|G (r,r )| 2 ~ — . (25) 



Lifting the ray density to the full phase space, Eq. (25) corresponds to rays moving out from 
the source in all directions according to 

cS(p-p ) , , 

Po{r,p;r ) = - r, 26) 

8noj\r — r Q \ 



where po = \p\(r — ro)/\r — r*o|. Note that from Eq. ([18]) , the ray dynamics take place on the 'energy' 
manifold c\p\ = 1. To obtain the source density on the boundary T, we set 

Po(r,p;r ) = pl(s,p s )S(c\p\ - 1), (27) 

where s parameterizes T and p s G denotes the momentum component tangential to T at s 

for fixed H(X) — c\p\ = 1; here -B^J 1 is an open ball in M.^ 1 of radius \p\ and centre s. Using 

local coordinates on the boundary, that is, p — (p s ,P ± ), then 8{p — po) = 6(p s — p Sn )b{p 1 ' — Pq) and 
|p|(5(p^ — Pq) — cp ± 5(c\p\ — 1). Combining these identities with (26) and (27) yields 



r, \ c 2 cos{0(r(s),r ))6(p s - p S0 ) 

p ( s »p s ; r o) = ^ i t n i , (28) 

OTTLj\r{s) — ro\ 

where 9(r(s),ro) is the angle between the normal to T at r(s) and the ray from ro to r(s). The 
resulting boundary layer density p^ is equivalent to a source density on the boundary producing 
the same ray field in the interior as the original source field after one reflection. 



The operator equivalent to the Frobenius-Perron operator (12), but now acting on boundary 
densities p T , is described in terms of a boundary integral operator B, 

Bp r (X s )= [ w(Y s ,u)S(X s - 1 (Y s ))p r (Y s )dY s (29) 

JdP 

where X s = (s,p s ), Y s = (s',p' s ) and 7 is the invertible boundary map. Also, <9P = Tx is the 

phase space on the boundary at fixed "energy" H(X) — c\p\. The weight w contains factors due 
to reflection/ transmission at interfaces within Q and a damping term of the form exp (—pL(s, s')), 
where L(s,s') is the length of the trajectory between s and s'. Note that, depending on the 
modelling assumptions, the factors in w may depend on uj. A similar weight term also needs to be 



applied in Eq. (28) since rays emanating from the source will also be subject to dissipation and 
possible reflection/ transmission at interfaces. We assume uniform damping in the interior, that 
is, the damping coefficient p does not depend on r. The damping term does then not alter the 
underlying ray dynamics |27j . Note that convexity is assumed to ensure 7 is well defined; non- 
convex regions can be handled by introducing a cut-off function in the shadow zone as in Ref. |21j 
or by subdividing the regions further. 

The stationary density on the boundary induced by the initial boundary distribution (X s , uS) 
can then be obtained using 

00 

/H = J2 B n {")pl{") = (i- #M)-Vo r M, (30) 

n=0 
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where B n contains trajectories undergoing n reflections at the boundary. 

The resulting density distr ibution on the boundary p r (X s ,ui) is then mapped back into the 
interior region. Using Eq. (23), one can relate p T to a wave energy density after projecting onto r 
space, that is for 57 C M 2 , 



p{r,u) = I p(r,p)dp, 

(31) 

p r (X*,u;)6(c\p\-l)\p\d\p\dQ, 



where O) is a polar coordinate system for p. One finally obtains 



p(r,w) = ± I p T (X s ,u;)dQ. (32) 



Now performing a change of variables from to s yields 

i \ 1 f r, Y s xCos(6»(r s ,r)) 

P\r,u) = -g / p {X ,w) — : : — ds, 33) 

c z J \r-r a \ 

where ?* s G T is the point on T in Cartesian coordinates corresponding to s. For systems with 
damping we must again apply a factor of the form exp (—pL(s, s')) to the final density evaluation 
(33). The total intensity \G(r, ro; lo)\ 2 may be estimated simply by adding \Go(r, ro; ui)\ 2 to the 
value of p(r, uS) given by ( [33] ) . 

We remark in passing that our approach is akin to the infinitesimal generator technique presented 
in [17j which aims at estimating the long term behaviour of the dynamics by numerically calculating 
the spectral properties of the generator of the Frobenius-Perron operator and extrapolating to long 
times. 



4. Numerical implementation: basis approximations and boundary element tech- 
niques 

The long term dynamics is contained in the operator (I — B)~ l and standard properties of 
Frobenius-Perron operators ensure that the sum over n in Eq. ( |30| ) converges for dissipative or 
open systems. In order to evaluate (I — *6) _1 , a finite dimensional approximation of the operator 
B must be constructed. In Refs. [23l [24] , basis expansions have been applied in both position 
and momentum coordinates, which is straightforward to implement for O C I 2 using univariate 
expansions in each argument. In particular, in Ref. [23] . a Fourier basis was employed globally on 
each subsystem in both position and momentum. It is well known that the convergence rate of a 
Fourier basis expansion depends on the smoothness of the approximant and also requires periodic 
boundary conditions |28j . The periodicity requirement will be fulfilled due to periodicity of the 
boundary curve and that the ray density corresponding to p s = ±|p| is zero giving periodicity in 
the momentum coordinate. However, the smoothness criterion often fails if, for example, corners 
are present in T; the initial density due to a point source is not then diffcrcntiable at the corners. 

In Ref. |24j it was suggested to split up the spatial approximation at corners in T and apply 
a basis approximation with no requirement for periodic boundary conditions in order to attain 
its optimal convergence properties. An approximation using a Chebyshev basis was implemented 
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since Chebyshev polynomials may be computed simply in terms of trigonometric functions and 
have similar convergence properties to Fourier basis approximations without the periodicity re- 
quirement. Such an approach also has a natural and optimal (in terms of polynomial degrees that 
can be integrated exactly) quadrature choice for the associated integrals, namely, Gauss- Chebyshev 
quadrature. The improvement in computational efficiency brought about by this approach opened 
up the method to larger scale multi-component problems. 

In this work we go further along this route by subdividing the boundary using a mesh. This 
introduces greater flexibility into the numerical method making it possible to both increase the 
polynomial order of the basis functions and to refine the mesh for the spatial coordinates. The 
advantage of refining the mesh is that additional degrees of freedom are introduced to the model 
whilst maintaining relatively low quadrature costs due to the lower order of the spatial basis ap- 
proximation. As such it is preferable to choose a basis which is orthogonal in the standard L 2 inner 
product. (From previous experience using a Chebyshev basis, where the associated orthogonal 
inner product is weighted, we found that low order approximations were very poor). Due to the 
mesh discretization procedure adopted here, the boundary conditions are not periodic and using a 
Fourier basis is not advisable; hence an alternative L 2 — orthogonal basis choice is proposed, namely, 
Legendre polynomials. 

The overall form of the proposed numerical scheme in 2D is therefore based on a Legendre 
polynomial approximation in both position and momentum augmented with a boundary mesh 
through which the spatial approximation may be localised. For an extension of the method to 3D 
cavities, see [29]. Recall that p s € (— \p\, \p\) and denote by p s a re-scaling of p s to (—1, 1). Let us 
also denote 

Pm(p.) = VV\p\P m {Ps)- 

An important special case of the scheme described above is when we fix the spatial basis approxima- 
tion at zero (constant) order and only refine the mesh, which essentially results in a scaled piecewise 
constant boundary element approximation. This type of approximation is also often referred to as 
Ulam's method |15U16j . although this would involve performing such an approximation in full phase 
space, rather than just in its spatial component. The implementation of separate meshes and basis 
approximations in both position and momentum would, however, complicate the procedure. In 
addition, our final computations will involve integration over spatial variables only. A spatial mesh 
therefore has the additional benefit of localising the integration regions, permitting the use of lower 
order quadrature rules. 

The overall approximation is then of the form 

n N s n p 

p r (X S ,L0) « E P*,l, m P a As)Pm(Ps), (34) 

a=l 1=0 m=0 

where N s and N p are the orders of the position and momentum basis expansions, respectively, and 
n is the number of elements in the mesh. Also 

P a ,l( s ) = VW^Pi(S a )x a {s), 

where s a parameterizes the a th element and is scaled to take values in (—1,1), \a is the char- 
acteristic function for element a and A a is the length of element a, a = l,..,ra. An analogous 
approximation is also made for p$ , the initial density, and the values of p a ,i, m are to be determined 
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by solving the resulting linear system. Note that the special case of piecewise constant bound- 
ary elements is reached by setting N s = above. The other limiting case, N p — 0, corresponds 
to a 'diffusive' reflection approximation, also denoted as Lambert's law [5], which is equivalent to 
the assumption that incoming rays are scattered uniformly in all directions after reflecting from a 
boundary. Equation ( 30 ) then becomes equivalent to the integral equations discussed in [STJ [3U] • 
Note that the case N p — N s = and n = 1 corresponds to an SEA approximation for a single 
subsystem (n equals the number of substructures for multi-component systems ) |23j . 



The matrix approximation B of the linear operator B ( 29 1 is given by taking the inner product 



< 



• > in 


L 2 (dP): 


B m +i- 


^N p (l+N 3 (a 


(2m 4 


l)(2i + l) 




4 


(2m 4 


1)(2J + 1) 




4 


(2m 4 


1)(2Z4 1) 



B(Pp,a(s')P b (p' s )) , P a ,l(s)Pm(Vs) 



w(Y s ,uj)5(X s - 1 (Y s ))P^ a (s')P b (p / s )P atl (s)P m ( Ps )dY s dX s 
w(Y\u J )P^ a (s , )P b (p , JP a Als(Yn)Prn(l P (Yn)dY s . 



(35) 



Here we write 7 — (7^,7^), to denote the splitting of the position and momentum parts of the 
boundary map. The boundary map can be very difficult to obtain in general, and hence we write 
the operator in terms of trajectories with fixed start and end points s' and s as follows 



B 



m+l+N p (l+N s {a~l)),b+l+N p (a+N s (l3-l)) 



(2m + l)(2Z + l) 



w(Y S : u)P m ( Ps (s : s'))P a As)Pb(p' s (s,s'))Pp,a(s') 



r 



dp' s 



ds 



ds' ds. 



(36) 



The resulting boundary integral formulation containing a pair of integrals over boundary coordi- 
nates bears a resemblance to standard variational Galerkin boundary integral formulations, see for 
example [3"T] . 



5. Numerical results 



In this section we consider a two-dimensional polygonal example whose boundary is meshed by 
subdividing each side into equidistant sections governed by a mesh parameter Ax. The number of 
elements on any given side is computed using the integer part of the side length divided by Ax. 



The Jacobian in Eq. ( 36 1 is written in the form 



dp' s 



\p\(n ■ (r — r'))(n' ■ (r — r')) 



-/|3 



(37) 



where n and n' are the inward unit normal vectors to T at r and r' , respectively. In addition r and 
r' are the Cartesian coordinate vectors corresponding to the boundary parameterised coordinates 



s and s' . In order to treat the corner singularities in (37), Gaussian quadrature is employed where 



end-points are not included as quadrature points. The convergence of the quadrature rules is still 
slow due to the peak in the integrand at corners. Telles' transformation techniques are employed 
to speed up the convergence [32]. 
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Figure 1: Five cavity polygonal example geometry 
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A five cavity system is considered as shown in Fig. [T] with Dirichlet boundary conditions on 
the outer boundary. The application of the method to multiple sub-domains is carried out using 
the same techniques as detailed in [24]. In general the configuration features irregular shaped, well 
separated pentagonal subsystems. In such a case one expects the wave energy to be fairly evenly 
distributed within subsystems and the energy levels in each subsystem to be distinct and well 
separated. The one feature of the configuration where such trends may not hold is due to the direct 
channel running to the right of the source point leading to a corridor type effect and energy being 
localised along this channel. Such dynamical features are likely to require higher order applications 
of our numerical algorithm to adequately resolve the behaviour. 

Energy distributions have been studied as a function of the frequency with a hysteretic damping 
factor r\ = 0.01, where fi = ujri/(2c) throughout tt. Here and in the remainder of this section the 
subsystems (or cavities) are numbered 1, 2, 5 from left to right. The wave speed is set to unity for 
simplicity. Note that the method allows c (and thus fj,) to vary between subsystems [24] , but we do 
not make use of this in the examples considered here. Therefore the subsystem interface reflection 



and transmission coefficients appearing in the weight term in (29) are simply and 1, respectively. 

Fig. [2] shows the energy distribution throughout the five subsystems compared with the results 
of a boundary element computation for the full wave problem. Explicitly we compute the square 
of the energy norm 

||G s: || 2 := / |G(r,r ;w)| 2 dr, £ = 1,2, ..,5. (38) 

The lines show the approximate energy norms given by the Liouville equation model computed for 
/ = (w/(27r)) = 10, 13, 16, 19,.. .,70 with N p = 10, N s = and n = 627, and hence the spatial 
approximation has been improved via mesh refinement rather than using higher order Legendre 
polynomial basis expansions. This has the benefit that the quadrature costs are kept moderate and 
leads to considerable efficiency gains. The computations here took approximately 2.5 hours per 
frequency compared with typical times of approximately 53 hours for a full tenth order Chebyshev 
basis computation using the methods reported in [24j . Here the computations were performed using 
a single core of a desktop PC with a 2.83 Ghz dual core processor. The circles each show the mean 
of 22 boundary element method computations, where the average is taken over a small range of 
frequencies centered around 10, 20,..., 70. The plus signs show the maxima and minima of the 
boundary element computations. The plot shows a good agreement between the solution of the 
Liouville equation and the mean of the full wave problem computations. In addition, the Liouville 
solution always lies within the range of the data for the full wave problem. 

Fig. [3] shows the convergence of our numerical algorithm in subsystem 5 as the precision of the 
approximation is increased by refining the mesh and simultaneously increasing the order of the 
Legendre basis approximation in direction space. The plot shows the relative error in the energy 
norm 1 1 1 1 2 compared against the value ||G§|| 2 computed with N p = 10, N s = and n = 627. The 
relative error plotted is thus given by ( 1 1 G\ \ \ 2 — 1 1 G5 1 1 2 ) / 1 1 G5 1 1 2 . The trend of the energy increasing as 
the precision of the numerical method is increased is shown by the error values all being positive and 
is expected since a low order implementation of our numerical scheme will not capture the direct 
energy channelled to the right of the source and into subsystem 5. Higher order basis functions 
in direction space are required to capture such directivity patterns in the wave field. In addition 
the plot shows the approximated solutions converging since the errors are getting closer to zero 
as higher order approximations are employed. For the case N p — 2, n — 70 one sees the error 
increasing as the frequency and thus damping are increased. For high damping, the direct channel 
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Figure 2: Energy distribution in five coupled acoustic cavities. Subsystem: 1 - blue (dash-dot), 2 - red (dashed), 3 ■ 
green (thin solid), 4 - black (dotted), 5 - cyan (thick solid) (colour online) 
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Figure 3: Evidence of convergence in subsystem 5 (colour online) 
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of energy (before reflection) from the source point will form a more significant part of the overall 
solution. For the case N p = 4, n = 70, the higher order basis in direction space has been able to 
model the channelling effect more accurately and the error has become almost uniform across the 
frequency range. 

Fig. [IJshows the energy density \G(r, r ;u>)\ 2 plotted for r G and / = 15. The upper subplot 
shows the energy density approximated using a boundary element method for the full wave problem. 
The central subplot shows the average of 22 full wave solutions at different frequencies 13 < / < 17. 
One sees that the averaging reduces the influence of small scale fluctuations in the solution and 
produces a generally smoother looking plot, showing more clearly the influence of the source point. 
The lower subplot shows the energy density approximated using the Liouvillc equation model with 
N p = 10, N s = and n = 627. Many similarities with the central subplot are apparent. For 
example, the direct energy channel to the right of the source and into subsystem 5, the location 
of the source and the dip in energy in the second subsystem close to the interface are features 
present in both plots. Also the energy channelled along the upper boundaries of subsystems 1 and 
4 is a common feature. This demonstrates that high frequency wave models based on the Liouville 
equation are indeed a powerful tool to obtain averaged solutions. Such an average solution can be 
obtained through averaging over an ensemble of similar structures or over a small frequency range. 
Applications are abundant in high frequency vibro-acoustics, where small changes in the domain 
due to, for example, manufacturing tolerances lead to large variations in the response and thus a 
simulation of any one particular problem will only be of statistical significance. Indeed, replacing 
highly oscillatory wave problems with smooth ray dynamical problems is often the only practical 
option in the high frequency regime. 

6. Conclusion 

The suitability of the Liouville equation as a model for high frequency wave problems has been 
discussed. A boundary integral approach has been proposed for its solution using the Frobenius- 
Perron operator. A numerical solution method has been implemented based on Legendre polynomial 
approximations in phase space together with a spatial boundary mesh. The proposed method has 
been verified by comparison with a BEM code for the full wave problem and numerical evidence of 
convergence has been shown. 
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Figure 4: Energy density compared against the full wave solution: upper subplot shows the full wave solution, central 
subplot shows the full wave solution averaged over a small frequency band, lower subplot shows the Liouville equation 
model (colour online) 



